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Abstract — The modulating functions method has been used 
for the identification of linear and nonlinear systems. In this 
paper, we generalize this method to the on-line identification 
of fractional order systems based on the Riemann-Liouville 
fractional derivatives. First, a new fractional integration by 
parts formula involving the fractional derivative of a modulating 
function is given. Then, we apply this formula to a fractional 
order system, for which the fractional derivatives of the input 
and the output can be transferred into the ones of the modulating 
functions. By choosing a set of modulating functions, a linear 
system of algebraic equations is obtained. Hence, the unknown 
parameters of a fractional order system can be estimated by 
solving a linear system. Using this method, we do not need any 
initial values which are usually unknown and not equal to zero. 
Also we do not need to estimate the fractional derivatives of noisy 
output. Moreover, it is shown that the proposed estimators are 
robust against high frequency sinusoidal noises and the ones due 
to a class of stochastic processes. Finally, the efficiency and the 
stability of the proposed method is confirmed by some numerical 
simulations. 

I. Introduction 

Fractional differential equations and fractional integrals are 
gaining importance in research community because of their 
capacity to accurately describe real world processes. The flow 
of fluid in a porous media, the conduction of heat in a semi- 
infinite slab, the voltage-current relation in a semi infinite 
transmission line are such examples of processes naturally 
modeled by fractional differential equations or fractional inte- 
grals. 

This paper is dealing with the identification of fractional 
order dynamical systems. The identification of such systems 
has been used for instance, for the estimation of the state of 
charge of lead acid batteries H], and for the identification in 
thermal systems ||2|, E). The goal of system identification is to 
estimate the parameters of a model from system input/output 
measurements. Different methods have been proposed for the 
identification of fractional order systems. Most of them consist 
in the generalization to fractional order systems of standard 
methods that were used in the identification of systems with 
integer order derivatives. We can classify these methods into 
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time domain methods and frequency domain methods. Time- 
domain methods have been introduced for example in ||4l, ||5], 
where a method based on the discretization of a fractional 
differential equation using Grunwald definition has been in- 
troduced, and the parameters have been estimated using least 
square approach. In |]6], a method based on the approxima- 
tion of a fractional integrator by a rational model has been 
proposed. In Q, the use of methods based on fractional 
orthogonal bases has been introduced. Other techniques can 
be also found for example in JS), ||9l, ifTOl and the references 
therein. 

In this paper, we are interested in the identification of 
fractional order systems using modulating functions method 
in case of noisy measurements. Modulating functions method 
has been developed by Shinbrot in ifTTI . lfT2l to estimate the 
parameters of a state space representation. Thanks to the prop- 
erties of modulating functions, the fractional differential equa- 
tion defining a fractional order system is transformed into a 
linear system of algebraic equations. Hence, instead of solving 
a fractional differential equation where the initial values are 
often unknown, the problem of identification is transformed 
into solving a linear system where the initial conditions are 
not required. Generalization of modulating functions method 
to fractional systems has been already proposed, for example 
in ifTsl . However, the authors of this paper proposed only to 
reduce the orders of the derivatives in a fractional differential 
equation. Moreover, the noisy case has not been considered. 

In the next section, basic definitions of fractional derivatives 
and modulating functions are recalled. Then in Section |lll] 
modulating functions method is applied to the identification 
of fractional order linear systems. Error analysis in both the 
continuous and discrete cases are presented in Section |IV] 
Numerical results are presented in Section |V] followed by 
conclusions summarizing the main results obtained. 

II. Preliminary 

In this section, first we recall the definitions and some useful 
properties of the Riemann-Liouville fractional derivative and 
modulating functions. Then, a new fractional integration by 
parts formula is given. 

A. Riemann-Liouville fractional derivative 

Let / be a continuous function defined on M, then the 
Riemann-Liouville fractional derivative of / is defined by (see 
lfT4l p. 62): Vi e K;, 
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where I — 1 < a < I with I G N*, and r(-) is the Gamma 
function (see ifTSl p. 255). As an example, using ([T]i the 
fractional derivative of an n*^ (n G N) order polynomial is 
given by (see HH p. 72): Vi G M+, 



(2) 



r(n + 1 - Of) 

We assume that the fractional derivative and the Laplace 
transform of / both exist, then the Laplace transform of the 
fractional derivative of / is given by (see |fT6l p. 284): V s e C, 

/-I 

£ {Dffit)} is) = sVis) [or'' VW] ' (3) 

1=0 

where / denotes the Laplace transform of /, and s denotes 
the variable in the frequency domain. 

Finally, we recall some results on the existence and the 
initial values of the fractional derivative in the following 
proposition. 

Proposition 1: (see HI PP- 75-77) If / G C'-i([0,T]) and 
G C{\0,T]) where T G M^, then the Riemann-Liouville 
fractional derivative D"/(t) exists, where I — \ < a < I with 
Z G N*. Moreover, the condition /('^O) = 0, for z = 0, . . . , 
1, is equivalent to the following condition: [D"./(t)]t^o ~ 

B. Modulating functions 

Let / G N*, T G M^, and g be a function satisfying the 
following properties: 

(Pi):.geC'([0,r]); 

{P2): 9^^K^) = 9^'\T) = Q, Vj=0,l,...,Z-l, 
then g is called modulating function of order I 1171 . Hence, 
according to Proposition [T] the modulating function g has also 
the following properties: 

(P3) : VO < /3 < T>1g{t) exists; 



(P4):V0 </?</, Df5(t) 



= 0. 



C. Fractional integration by parts 

Another useful property of the modulating functions is given 
in the following theorem. 

Theorem 1: Let y be a function such that the a*'' order 
fractional derivative exits and g be a modulating function of 
order I with I — 1 < a < I with / G N*. Then, we have: 

/ g{T^t)'D'iy{t)dt= [ D'^ g{t) y{T ~ t) dt, (4) 
JO Jo 

where T eR*^_. 

Proof: By applying the convolution theorem of the 

Laplace transform (see lITSl . p. 1020), we get: 



Then, using (O we obtain: 

g{s)C{Bfy{t)}{s)^ 
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Moreover, using (|3) and (P4) we get: 

£{D'^g{t)}is)=s'^g{s), (7) 

C{g^^\t)}{s)^s^g{s), (8) 

for i ^ 0, . . . ,1 — 1. Consequently, by applying (|7]i, (|8]l and 
the inverse of the Laplace transform to (|6), we obtain: 

C-'{gis)C{D?y{t)}is)}iT)^ 

C-'{C{B?g{t)}{s) y{s)}{T) 



(9) 
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Using {P2), the initial conditions [D"~'^~^y{t)]^_^, for i = 
0,1, ... ,1 ~ 1, can be eliminated. 

Finally, this proof can be completed by applying the con- 
volution theorem of the Laplace transform to (|9). ■ 

According to the previous theorem, we can see that by 
working in the frequency domain: on the one hand, we can 
obtain an integral formula which can be considered as the 
generalization of the classical integration by parts formula. 
Another fractional integration by parts formula has also been 
given in ITSl , however the Caputo fractional derivative was 
involved in the formula; on the other hand, the initial con- 
ditions of the fractional derivatives of y can be eliminated 
using a modulating function. In fact, the idea of obtaining 
this theorem is inspired by the recent algebraic parametric 
estimation technique JTU, JSO), lEI], El, ||23l, ||24l, IBll . 
Il26l . which eliminates the unknown initial conditions by 
applying algebraic manipulations in the frequency domain. 

III. Identification of fractional order systems 

A. Fractional order linear systems 

In this section, we consider a class of fractional order 
linear systems which are defined by the following fractional 
differential equation: Vt G / = [0, T] C R^^ 



M 

3=0 



(10) 



where y is the output, u is the input, a,; , bj G R* are unknown 
parameters to be identified, and , /3j G R+ are assumed 
< ao < ai < ■ • • < aL, < /3o < /?! < ■ • • < ;Sj\f with 
L, M G N. 

Let y^ be a noisy observation of y on the interval /: 

ytel, y^{t) = y{t)+mit), (11) 

where zu is an integrable nois^H We are going to estimate 
the unknown parameters in (fTol i using the input u and the 
observation y^. 

One of the standard methods for the identification of 
systems with integer order derivatives is to use the least- 
squares method ll28l . This method was generalized to frac- 
tional order systems in ||4l, ||5]. In this method, we need 

'More generally, the noise is a stochastic process, which is integrable in 
the sense of convergence in mean square 1221 . 



to estimate the fractional derivatives of y and u using a 
fractional order differentiator [E9l . [l30l . llSTI . However, a 
fractional order differentiator often contains a truncated term 
error An alternative method consists in solving the problem 
in the frequency domain by applying the Laplace transform to 
( [Tol l. However, according to (|3]l, this application can produce 
unknown initial conditions. In the next subsection, we are 
going to apply modulating functions method to eliminate these 
unknown initial conditions. 

B. Application of modulating functions method 

We denote W = L + M + 1, a = max (ai, /3a/), and 
I = \a\ , where \a\ denotes the smallest integer greater than 
or equal to a. Then, we take a set of modulating functions 
{5n}^=i ^i'^h < € N. Using Theorem [T] we can give 
the following proposition. 

Proposition 2: Let {gn}^^i be a set of modulating func- 
tions of order I. If we assume that 6o = 1 in the fractional order 
linear system defined by ( fTOl i. then the unknown parameters in 
this fractional order linear system can be estimated by solving 
the following linear system: 
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where B ^ \^bi, ■ ■ ■ , bujj , A = {gq, ■ ■ ■ , gl) are estima- 
tors of the unknown parameters, and 

UN{n,j) = - f B^^^g^{t)u{T~t)dt, (13) 
Jo 



[ BT9n{t)y^{T~t)dt, (14) 



lN{n) = I B^°gn{t)u{T-t)dt, 
'o 



(15) 



for n = 1, . . . , iV, j 1, . . . , A/ and i = 0, . . . , L. 

Proof: By multiplying the modulating functions 5„ to the 
equation ( flOl i and by integrating between and T, we get: 



rT L 



•^0 z=0 

T M 

/ y2bjgn{T-t)D^^u{t)dt, 



(16) 



are known and without noise. On the other hand, if the 
fractional derivatives of u cannot be analytically calculated 
or are difficult to calculate, we can solve the problem by 
calculating the ones of the modulating functions. 

Finally, let us mention that since the proposed estimators are 
given in causal case, if we take the value of T to be equal to the 
time where we estimate the parameters, then these estimators 
can be used for on-line identification applications. 

IV. Error analysis 

In this section, we are going to study the noise effect in the 
integrals obtained in Proposition |2] For this purpose, we study 
the noise error contributions due to a high frequency sinusoidal 
noise and the ones due to a class of stochastic processes in 
continuous case and in discrete case, respectively. 

A. Error analysis in continuous case 

There are many applications where the output signal is 
corrupted by a sinusoidal noise of higher frequency ll26l . 
Hence, we assume that the noise w is a high frequency 
sinusoidal noise in this subsection. 

By writing y^ ~ y + vd, the integral {n,i + l) given in 
(O, for n ~ 1, . . . , N and i ^ 0, . . . , L, can be divided into: 

Y^in,i + l)=YN{n,i + l) + e^in,i + l), (18) 

where Y/v(ri, i + 1) is given in ( [TtI i. and the associated noise 
error contribution is given by: 

{n,i + l)= [ D^^gnit)nj{T -t)dt. (19) 
Jo 

Consequently, the estimation errors for the estimators given in 
Proposition [2] only come from these noise error contributions. 
In the following proposition, we give error bounds for these 
noise error contributions. 

Proposition 3: We assume that Vt G I, w{t) ~ csm{ujt + 
(p) with c, w G M.*^_ and (f> £ [0, 27r[. Moreover, we assume 
that D"*^^g„(<) exists and is continuous on [0,T], for n = 
1, . . . ,N and i — 0, . . . , L. Then, we have: 

\e^in,i + l)\ < -ra.+i + -|[D^'g„(t)],^j,|, (20) 

id LO 

where Cq;+i = sup |D"'+^(7„(t)|, and e^(ri-, i + is given 



te[o,T] 



for n = 1, . . . , A^. Then, using Theorem [T] we get 

M L „T 

j = l 1=0 J^ 

(17) 

where J7jv(n, j) = - /(f of^ c/„(t) u{T - t) dt, and 
YN{n,i + l) = /(f B^'gn{t) y{T - t) dt. Finally, this proof 
can be completed by substituting y by y^ in (fTTI i. ■ 
Consequently, thanks to Theorem [T] instead of estimating 
(resp. calculating) the fractional derivatives of y (resp. it), we 
calculate the fractional derivatives of the modulating functions. 
On the one hand, comparing to y, the modulating functions 



by dUll. 

Proof: Since ^ {D°"g„(i)} = D°'+^g„(i) exists, then 
by applying integration by parts and (P4), we get: 



e^(n,i + l)= / D'^'gn{t)csui{uj{T -t) + (l))dt 
Jo 

D""+\g„(t) cos(a;(T - t) + 0) dt 

lo 



(21) 



If D^'+^g„(i) e C([0,T]), then this proof can be completed 
using ( l2n i. ■ 



According to the previous proposition, if the frequency of 
the sinusoidal noise is high, then the associated noise error 
contributions can be negHgible. Consequently, the estimators 
given in Proposition |2] can cope with this kind of noises. 

B. Error analysis in discrete case 

From now on, we assume that the noisy observation 
defined in (fTTl i is given in a discrete case. Let y^{tj) = 
y{tj) +w{tj) be a noisy discrete observation of y given with 
an equidistant sampling period Tg, where = £ f^* ^ 

and tj = jTs, for j = 0, - ■ ■ , m. 

Since y'^ is a discrete measurement, we apply a numerical 
integration method to approximate the integrals in (fTZt . Let 
wq > 0, Wm > and Wj > for j = 1 , . . . , to — 1 be 
the weights for a given numerical integration method, where 
weight qq (resp. is set to zero when there is an infinite 
value at <o = (resp. tm = T). Then, the integral Y^{n^ i + 
given in ( fT4l i. for n — 1, . . . ,N and i = 0, . . . , L, can be 
approximated by; 

m 

y^„(n, z + 1) := T, ^ w, .9("*)(t,) r (i™-,), (22) 
i=o 

where gl^'\tj) := [^t' Sn{t)]t=t -- The integrals given in (fTST i 
and ( fTSl l can be approximated in a similar way. 
By writing y^{tj) = y{tj) + w{tj), we get: 

Y^^„^{n,i + 1) = YNMri,i + I) + e^^,^{n,i + I), (23) 

where 

m 

YN,m{n, i + l)=TsY, ^1 5i"-^ {tj) y{tm-j), (24) 

j=0 
m 

e^^^in, * + 1) = T, ^ w, gtHt,) w{t,^-,). (25) 

Thus the integral Y^{n, i + 1) is corrupted by two sources of 
errors: 

• the numerical error which comes from a numerical inte- 
gration method, 

• the noise error contribution i + 1). 
Consequently, the estimation errors for the estimators given 
in Proposition [2] come from both the numerical errors and the 
noise error contributions in the discrete noisy case. 

It is well known that if the value of T is set, then when Tg 
tends to 0, i.e. m — > +oo, the numerical errors tend to 0. In 
the next subsection, we are going to study the effect of the 
sampling period on the noise error contributions. 

C. Influence of sampling period on noise error contributions 

In this subsection, we consider a family of noises which are 
stochastic processes satisfying the following conditions: 
(Ci) : for any t,s G I, t ^ s, w{t) and w{s) are 
independent; 

(C2) : the mean value function of w{-) denoted by E[-] 
belongs to C{I); 



(C3) : the variance function of w{-) denoted by Var[-] is 
bounded on /. 

Note that the white Gaussian noise and the Poisson noise 
satisfy these conditions. Then, we can give the following 
proposition. 

Proposition 4: Let ti7(-) be a stochastic process satisfying 
conditions (Ci) — {C3), and vuitj), for j = 0, • • • , m, be a 
sequence of {w{ )} with an equidistant sampling period Tg. 
If D"'5„(f;) g then we have the following convergence 

in mean square of the noise error contribution in the integral 
Y^{n, i + given in (O, for n = 1, . . . , and i = 0, . . . , L: 

e^.„(n, I + 1) £^ j ^Tgnit) E [w{T - t)] dt, (26) 

where e^,„(ri, i + 1) is given in ( |25] |. Moreover, if Vt G /, 
E [vj{t)] = 0, then we have: 

e^^„K* + l)=^0. (27) 

The proof of the previous proposition can be obtained in 
a similar way to the one given in [i27l. Moreover, a similar 
result was studied using the non-standard analysis in ll32l . 

Consequently, according to the previous proposition, the 
noise error contributions can be increasing with respect to the 
sampling period. 

Finally, let us mention that solving the linear system given 
in Proposition |2] in noisy case is related to the matrix pertur- 
bation theory. The accuracy and the stability of the proposed 
estimators not only depend on the noise error contributions, 
but also depend on the condition number of the associated 
matrix. This condition number depends both on the input 
and the output of the fractional order system and on the 
used modulating functions. In general, we should choose the 
modulating functions that can give a small condition number. 
This study is out of the scope of this paper. 

V. Simulation results 

In order to illustrate the accuracy and robustness with 
respect to corrupting noises of the proposed estimators, we 
present some numerical results in this section. 

Let us consider a fractional order system defined by the 
following fractional differential equation: \/t G [0,8], 

aoB^°y{t) + aiD^Hj(t) + a^B'^^yit) = u{t), (28) 

where cq = 3, ai = 2, 02 — 1, aQ — 0, ai ~ 0.5 and 
a2 = 1-5. We assume that the output is y{t) = sin(3t) + 1. 
Hence, the initial conditions of y are not equal to 0. Moreover, 
the expression of the input can be obtained using (|2]i and the 
following formula (see [[331 p. 83): 

D^' sin(3t) = 

J!^Z_ ( 1; i(2 - a,), -{3 - a,); --ShA , ^^'^^ 
r(2 - ai) V 2^ ' 2^ ''4 y ' 

where pFq {ci, . . . ,Cp;di, . . . ,dq; •) is the generaHzed hyper- 
geometric function given in 133] p. 303. 



In our identification procedure, we use the following mod- 
ulating functions, the fractional derivative of which are simple 
to calculate: 

where n = 1,2, ■ ■ ■ , N with N = 13. The value of T is taken 
to be equal to the time where we estimate the parameters. Let 
us recall that this kind of functions has been obtained when 
the algebraic parametric estimation technique was applied to 
the parameter estimation for signals described by differential 
equations II2TI . 

In the two following examples, we estimate the parameters 
ao, fli and 02 using the noisy observation of y where the 
noise is a high frequency sinusoidal noise and a gaussian 
noise, respectively. Moreover, we apply the trapezoidal rule 
to numerically approximate the integrals in our estimators. 

Example 1. In this example, we assume that y^{tj) = 
y{tj)+0.5 sin{10Hj) with = 0.01. We can see this discrete 
noisy signal in Figure [T] In our identification procedure, we 
take T = ti for ti G [1.8,8]. The obtained estimations and 
the associated relative estimation errors are shown in Figure |2] 
and Figure [3] Hence, we can see that the proposed estimators 
can cope with a high frequency sinusoidal noise. 

Example 2. In this example, we assume that y^{tj) ~ 
y{tj) + <jw{tj), where Tg = 0.01, ■Dj{tj) is simulated from 
a zero-mean white Gaussian iid sequence, and a G is 
adjusted in such a way that the signal-to-noise ratio is equal 
to SNR = 22dB. This noisy observation is shown in Figure|4] 
The obtained estimations and the associated relative estimation 
errors are given in Figure 2] and Figure We can see that the 
proposed estimators are robust against a gaussian noise. 




Fig. 1. The noise-free output and its noisy observation. 

VI. Conclusions 

In this paper, the modulating functions method has been 
generalized to the on-line identification problem of fractional 
order systems. Using this method, the unknown parameters 
have been estimated by solving a linear system of algebraic 
equations involving the input and the noisy output. Thanks 
to the properties of modulating functions, we do not need to 
estimate the fractional derivatives of the output, to calculate 
the ones of the input. We do not need to know the initial 
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Fig. 2. The exact parameters and their estimations with a sinusoidal noise. 
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Fig. 3. Relative estimation errors in the sinusoidal noise case. 
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Fig. 4. Tlie noise-free output and its noisy observation. 



conditions either Moreover, the integral given in the estimators 
can reduce the noise effects due to a high frequency sinusoidal 
noise or a class of stochastic processes. The efficiency and the 
robustness against corrupting noises have been confirmed by 
numerical examples. In order to improve the robustness against 
noises, some methods such as the instrumental variable method 
will be applied ll34l . It was mentioned that the proposed 
estimators also depend on the choice of modulating functions. 
This problem will be studied in the future work. Moreover, the 
estimation of the fractional derivative orders in a time-delayed 
fractional order system will be considered. 
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Fig. 5. The exact parameters and their estimations with a gaussian noise. 




Fig. 6. Relative estimation errors in the gaussian noise case. 
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